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Abstract 

We study S'[/(3)-breaking effects in the neutral B^-Bd and Bg-Bg systems with unquenched 
^ ■ Nf = 2 + 1 lattice QCD. We calculate the relevant matrix elements on the MILC collaboration's 

d . gauge configurations with asqtad-improved staggered sea quarks. For the valence light-quarks {u, 

d, and s) we use the asqtad action, while for b quarks we use the Fermilab action. We obtain 
^ = fBs\/BBs/fBci\/BBi = 1-268 ± 0.063. We also present results for the ratio of bag parameters 
BbJBb^ and the ratio of CKM matrix elements |Vtrf|/|Vts|- Although we focus on the calculation 
of ^, the strategy and techniques described here will be employed in future extended studies of the 
B mixing parameters AM^^^ and AF^j^s in the Standard Model and beyond. 
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I. INTRODUCTION 



The observation of new particles at high-energy colhders is not the only way for new 
physics to be discovered. It can also be unveiled through the observation of deviations 
from the Standard Model (SM) via high-precision measurements of low-energy observables 
in high-luminosity experiments. This requires matching precision in the theoretical SM 
predictions for these observables. In principle, such a comparison could reveal the exchange 
of virtual, new heavy particles involving scales much higher than those that can be achieved 
in direct production at high-energy colliders. 

Heavy-flavor physics and, in particular, neutral-meson mixing are potentially very sen- 
sitive to these virtual effects. Neutral-meson mixing occurs at loop level in the SM, see 
Fig. [H and it is further suppressed by small Cabibbo-Kobayashi-Maskawa (CKM) matrix 
elements, so the effect of new particles in the internal loops could be noticeable in the pa- 
rameters describing the mixing. Indeed, there are several measurements for which there is 
a 2-3cr difference from the SM prediction. These include sin(2/3) the like-sign dimuon 
charge asymmetry [2|, and unitarity triangle (UT) fits j3|-0]. It has been argued that these 
differences may be due to physics beyond the Standard Model (BSM) affecting the neutral 
i?- meson mixing processes [3|, IJ] . 

In the Bg system, the relative phase between the decay amplitudes with and without 
mixing, /J^, could also show BSM effects, as pointed out in Ref. jl] and later hinted at in 



a Tevatron measurement [9|. Although new measurements at CDF [13] and D0 11 are 
in better agreement with the SM, reducing the difference from ~ Sa to ~ la, there is still 
room for a large deviation of Ps from SM values. 

The main parameters describing mixing in the and the B^ systems are the mass 
differences, AMs(rf), and the decay width differences, AFg^^i), between the heavy and light 
mass eigenstates, and the CP violating phases (psid)- The phases (j)s{d) are defined as 
the argument of the ratio of the dispersive and absorptive off-diagonal elements of the time 
evolution matrix which describes the mixing [l^ . The existence of new, heavy particles in 
loops could affect the value of the mass differences, given by the dispersive part of the time 



evolution matrix. The mass differences AMg jl3Hl5| and AM^ [16] have been measured 
with an accuracy better than 1%. Improving the theoretical control on these quantities is 
thus crucial in order to fully exploit the potential of CP violating observables to search 
for nonstandard physics. In addition, the theoretical calculation of BSM contributions to 
mixing and the experimental measurement of -B° mixing parameters can help in constraining 



BSM parameters and understanding new physics [6|. Several recent studies have addressed 



that task [3|-[7|, |17H24||. finding that one of the main limitations to further constraining the 
parameter space in BSM theories is the error associated with the theoretical calculation of 
the nonperturbative inputs. 

The most interesting quantity to analyze in B^ mixing phenomena is the SU (3)-breaking 
ratio ^, which measures the difference between the mixing parameters in the B^ and the 
B^ systems, and enters the relation between the ratio of mass differences and CKM matrix 
elements as 



Vt 



td 



Vt 



ts 



Its value, together with the experimental measurement of the mass differences AM^ deter- 
mines the ratio of CKM matrix elements | Vt^/Visl, which constrains one side of the unitarity 
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FIG. 1. Box diagrams contributing to — mixing in the SM. Gluon exchanges shown in the 
plot are just representative of the QCD corrections. 



triangle 25, 26|. Thus, ^ is one of the key ingredients in UT analyses js], M 

In the SM, mixing is due to box diagrams with the exchange of two H^-bosons, like those in 
Fig. [H These box diagrams can be rewritten in terms of an effective Hamiltonian with four- 
fermion operators describing processes with AB = 2. In BSM theories, mixing processes can 
receive contributions from additional diagrams due to the exchange of new, heavy particles. 
These can also be parametrized in terms of four-fermion effective operators built with SM 
degrees of freedom. The most general effective Hamiltonian describing processes with AB = 
2 was given m and can also be found in [29|. There are a total of five independent 



operators (plus parity conjugates) in the Hamiltonian, but only three of them contribute to 
mixing in the SM 



"^rfF^M — ^ CiOi , (1.2) 



1=1 



with 



Ol= {q^LV) [q^LU) , (1.3) 

where i and j are color indices, and L and R are the Dirac projection operators |(1 — 75) 
and ^(1 +75) respectively. The fields q denote strange or down fields for and B^ mixing 
respectively, and h represents the bottom field. 

The matrix element of the first operator in Eq. (11.31) . 0\, provides the mass difference in 
the SM: 



AMf^ = ^^jr^\V:M"^2So{xt)MsjlBs, , (1.4) 



67r2 

where S'o(xt) is the Inami-Lim function [30|, which depends on the top quark mass through 
Xt = m^/M^, and the quantity r^^ is a perturbative QCD correction factor. The products 
f^^BBg parametrize the hadronic matrix elements in the effective theory by 

{B',\Ol\B'^){f,) = lMljlBs,{f^). (1.5) 

The factors /^^ are the B^ decay constants. The renormalization group invariant bag pa- 
rameters in Eq. (11.41) are related to the scheme and scale dependent bag parameters in 
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(11. 5p at next-to-leading order (NLO) by 
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where J5 is known in both MS-NDR (naive dimensional regularization) and MS-HV ('t Hooft- 



Veltman) schemes 3l|]. Bag parameters have traditionally been used to measure the de- 



viation of the four-fermion operator matrix elements from their vacuum insertion values, 
Bb = 1. 

The S'f/(3)-breaking parameter ^ can be written in terms of decay constants and bag 
parameters 



fss \/ Bbs 

fsa a/ Bb^ 



(1.7) 



Many of the uncertainties that affect the theoretical calculation of the decay constants and 
bag parameters cancel totally or partially in this ratio, leaving the chiral extrapolation as 
the dominant error. Hence, ^ and the combination of CKM matrix elements related to it, 
can be determined with a significantly smaller error than the individual matrix elements. 

The hadronic matrix elements in Eq. (II. 5p encode the nonperturbative physics of the 
problem and are best calculated using lattice QCD. Our current knowledge of them limits 
the accuracy with which the CKM matrix elements appearing in Eq. (II. 4p can be determined 
from the experimental measurements of AMg^d)- In particular, the uncertainty associated 
with the calculation of ^ is one of the main limiting factors in UT analyses, so improvement 
in the knowledge of ^ is crucial to disentangle the origin of the 2-3o" tension. 

There are two 2 + 1 unquenched lattice calculations of the ratio in the literature. One 
is by the HPQCD collaboration [32|, which quotes the value C = 1.258(33). The other is an 
exploratory study by the RBC and UKQCD collaborations [33| on a single lattice spacing 
and using the static limit for the bottom quark; their result is ^ = 1.13(12). In this paper, 
we report a lattice calculation of ^ at the few percent level. 

In Ref . 



Preliminary results related to the work here were presented in |34l437 



34 , the 



simulation and correlator fitting methods were described using data for one lattice spacing, 
while Refs. [35|, 36 1 focused on the discussion of statistical and fitting errors, and the chiral 
extrapolation method. In Ref. 37| we studied the matching method and the heavy-quark 



discretization errors. 

The primary difference between this work and the HPQCD calculation in Ref. [32 
the treatment of the valence b quarks. The HPQCD collaboration uses lattice NRQCD 
while we employ the clover action 



with the Fermilab interpretation |40 



An advantage 

of the Fermilab method is that it can also be efficiently used to simulate charm quarks, so 
the analysis performed in this work can be easily extended to the study of the short-distance 
contributions to D^-D^ mixing. Although in the case of neutral D mixing the long-distance 
contributions are believed to be dominant, a calculation of the short- distance contributions 
nevertheless can provide valuable constraints on extensions of the SM (4l| . 

In order to achieve the few-percent level of precision required by phenomenology, we use 
lattice QCD simulations with realistic sea quarks. In particular, ^we employ a subset of the 
MILC configurations with 2+1 fiavors of asqtad sea quarks liT 



44 



In the valence sector, 

we use the same staggered asqtad action to simulate the hght quarks. The configurations 
we use in this analysis were generated using the fourth-root procedure for eliminating extra 
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degrees of freedom originating from fermion doubling. Despite the nonlocal violations of 

there are strong theoretical 
that the local, 



unitarity of the rooted theory at non-zero lattice spacing [45|, |46 



arguments |47H50l] . as well as other analytical and numerical evidence |51l-l54 
unitary theory of QCD is recovered in the continuum limit. This gives us confidence that 
the rooting procedure yields valid results. We also explicitly tested the rooting procedure 
as well as improvements in our heavy action by calculating the spin-dependent hyperfine 
splittings for Bg and Dg mesons in Ref. jsHj. 

Our collaboration has already successfully used the asqtad MILC ensembles in similar 
calculations of other quantities involving B mesons, as part of a broad program of calculating 
matrix elements: for example, the extraction of the CKM matrix elements \Vub\ and \Vcb\ 
from the calculation of, respectively, the semileptonic form factors describing the processes 



56| and B — >■ 0*1^ [57|, |58|; or, more recently, the calculation of the fs and /^^ 



59| and the form-factor ratios between the semileptonic decays B — j- D~^l v 



60 



B -nlv 
decay constants 
and Bs ^ 

This paper is organized as follows. In Sec. [ITl we describe the actions and parameters 
used in our numerical simulations, as well as the construction of the mixing operators and 
correlation functions. Section UTT] presents the renormalization method using one-loop mean- 
field improved lattice perturbation theory. We include a discussion of the errors associated 
with the matching and numerical values of the matching coefficients used. Next, in Sec. IIV[ 
we give the details of the procedure for the correlator fits. Section |V] is devoted to the chiral- 
continuum extrapolation, which is performed within the framework of rooted staggered chiral 
perturbation theory (63-67 



We describe and discuss the choice of the functional form used 
in the extrapolation, the different fitting methods tested, and the choice of parameters and 
parametrization. In Sec. IVIj we list and estimate the different systematic errors. Finally, 



Section [VIII compiles our final results for the parameter ^ as well as for |Vjrf|/|Vts|; and the 
ratio of bag parameters BbJBb^- We also discuss planned future improvements in the 
calculation of B^ mixing parameters by our collaboration. In Appendix |Al we provide the 
explicit formulas for the chiral fit functions used in the chiral fits described in Section |Vl 
In Appendix |Bl we compile the functions needed to estimate the heavy-quark discretization 
errors in our calculation. Finally, Appendix O discusses our choices for prior central values 
and widths for the correlator fits. 



II. NUMERICAL SIMULATIONS 
A. Parameters of the simulations 

The n/ = 2 + 1 MILC ensembles [el used in our calculation include the effect of three 
sea-quark flavors: two degenerate light quarks corresponding to the up and down quarks 
(although with larger masses than the physical ones), and one heavier quark corresponding 
to the strange quark. These dynamical quarks are simulated using the asqtad improved 
staggered action with errors starting at 0{asa^) j68|. The gluon action is a Symanzik 
improved and tadpole improved action, with 0{asO?) errors coming from the gluon loops 
removed jsH, The couplings needed to remove the 0{asO?) errors coming from quark 



loops [71[ were available only after the generation of conflgurations was well advanced, so 
these effects are not accounted for in the MILC ensembles. The dominant errors in the 
gauge action are thus also of 0{a^, agO?)- 

The valence light-quark propagators are generated using the asqtad action and converted 



5 



TABLE I. Parameters of the ensembles analyzed in this work. The first two rows show the approx- 
imate lattice spacing and the volume, ami and anih are the light and strange sea quark masses, 
respectively. A'confs is the number of configurations analyzed from each ensemble, and aniq are the 
light valence quark masses. The ri/a values are obtained by fitting the calculated ri/a to a smooth 
function [6l|, as explained in Ref. [g^ ]. 



^ a(fm) 


(-)'x - 


anil / arrih 


-^confs 




arriq 


ri/a 


0.12 


243 X 64 


0.005/0.05 


529 


0.005, 


0.007, 0.01, 0.02, 0.03, 0.0415 


2.64 


0.12 


20^ X 64 


0.007/0.05 


833 


0.005, 


0.007, 0.01, 0.02, 0.03, 0.0415 


2.63 


0.12 


20^ X 64 


0.01/0.05 


592 


0.005, 


0.007, 0.01, 0.02, 0.03, 0.0415 


2.62 


0.12 


20^ X 64 


0.02/0.05 


460 


0.005, 


0.007, 0.01, 0.02, 0.03, 0.0415 


2.65 



0.09 28^ X 96 0.0062/0.031 557 0.0031, 0.0044, 0.062, 0.0124, 0.0272, 0.031 3.70 
0.09 28^ X 96 0.0124/0.031 534 0.0031, 0.0042, 0.062, 0.0124, 0.0272, 0.031 3.72 



to naive quark propagators using the relation [72 

eered 



(2.1) 



Xq X\ X2 xs 

7o 7i 72 73 • 



where Q{x) 

For the heavy bottom quarks we use the Sheikholeslami-Wohlert action 



39 with the 



Fermilab interpretation for heavy-quark systems 4^]. This interpretation retains the full 
mass dependence of the theory within the parameters of the lattice action. A tree-level 
matching to QCD is then performed via heavy quark effective theory (HQET), after which 
it can be shown that the errors in the action begin at OlagAQCBd, -^qcd^^"^) times bounded 
functions of m^a, the 6-quark mass in lattice units. 

We perform our analysis at two different values of the lattice spacing, a ~ 0.12, 0.09 fm, 
and for a variety of sea-quark masses. The values used are shown in Table [H The mass of 
the heavy b quark is fixed to its physical value by computing the spin-averaged Bg kinetic 
This determines the b quark's hopping parameters 
0.0923 for the 



mass 



55 



Kb 



0.09 fm lattice 55 



- 0.0860 for the a ^ 
and thus the bare b 



0.12 fm lattice and 

quark mass. We simulate the B mesons with the six different values of light- valence quark 
mass listed in Table [H the smallest of which is around rris/S, in order to facilitate the 
extrapolation/interpolation to the physical down/strange quark masses. 



B. Correlators: the open-meson propagator 

As described in the introduction, the study of the S't/(3)-breaking ratio ^ requires the 
calculation of the hadronic matrix element^ (C^i) and (Of), the latter of which mixes with 
{Of) under renormalization, for both q = d,s. The matrix elements are obtained from 
three-point correlation functions with zero spatial momentum 

Co^{t^,ty) = J2{B'q{ty,y)OmB'q{tx,xy) , (2.2) 

x,y 



1 To simplify the notation, wc define (Of) = {B'^\0I\B'^) for i = 1,2. 
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where i = 1 or 2, and the S-meson creation operator B^(t, xY = b(t, x')S{x, x')'-f5q{t, x), 
with q{t, x) the naive hght quark field, whose propagator is constructed from the staggered 
propagator via Eq. (12. ip . and with S{x,x') a smearing function. Our choice of smearing 
function is discussed in Sec. IIV A[ The structure of the functions in f l2.2p is depicted in Fig. |2J 
The four-fermion operators are placed at the origin while 5-mesons are positioned at 
X and y. This layout allows us to perform the three-point function fits over both t^. and 
ty, maximizing the information included in the fits. In order to extract the relevant matrix 
elements from (12. 2p . we need to determine the overlap of the 5- meson creation operator 
with the ground state. Therefore, we also need the pseudoscalar two-point correlator with 
zero spatial momentum 

Cp5(t) = $^(5;(t,a;)5j(0,0)t). (2.3) 

X 

The calculation of both three-point and two-point correlators can be organized into con- 
venient structures. Starting with a general correlator with Dirac structure Fi x F2, which 
accommodates a full set of AB = 2 four-quark operators, including those in Eq. (II. 3p . 

Csit., ty) = $^(5j(t„ 2/)g(0)Fi6(0)g(0)F26(0)Sj(t., x)^) , (2.4) 

x,y 

and performing the four possible Wick contractions, we obtain 

C3{t,,ty) = 5^{tr[75L,(x,0)Fii7,(0,x)]tr[75L,(y,0)F2i/b(0,y)] 

+ tr[75L,(y, 0)T,H,{0, y)] ti[^,Lq{x, 0)TMO, x)] 

- tr[75L,(x, 0)Fii7fe(0, yh5L,{y, 0)F2i/b(0, x)] 

- tr [75L,(x, 0)F2i/fe(0, yHL.iy, 0)T,H,{0, x)] } (2.5) 
= J2{ 0)Fi75i/J(x, 0)] tr[L,(y, 0)F275^^;(y, 0)] 

x,y 

+ tr[L,(y, 0)Ta5Hliy, 0)] tr[L,(x, 0)F275i^J(x, 0)] (2.6) 

- tr[L,(x, 0)Ta5Hl{y, 0)L,{y, 0)T2J,hI{x, 0)] 

- tr[L,(x, 0)T2i5Hl{y, 0)L,{y, 0)T,^,hI{x, 0)] }, (2.7) 

where Lg is the (naive) light-quark propagator, and Hh is the heavy-quark propagator. The 
traces in Eq. (12.51) run over spin and color indices. 

These correlators can be rewritten as 

Cs{t., ty) = Ff (t.)Fri?r;(t,) + Ff (t,)Fri?r;(tx) (2.8) 
-Ff i^r (^.)rr^ca^(i.) - Ff i^r (^.)rri^ca^(^.) , 

where summation over repeated indices is implied, and we have introduced the basic objects 

K^t) = irKfait, 0)L;i(t, 0) , (2.9) 

with Dirac indices labeled as a,/3,a,T and color indices labeled as a,c,d. We call the 
combination of propagators £'^^(ta;) defined in Eq. (12. 9p "open-meson propagator" . Once the 
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to (source) 



FIG. 2. Structure of the three-point correlators. A is created at rest at tx + to < to- At time 
to, it osciUates into a via the operator O?, which is subsequently annihilated at ty + to > t^. 



open-meson propagators have been computed and saved, all correlation functions needed for 
i?-meson mixing, including BSM operators, can be immediately constructed by contracting 
them with the appropriate Dirac structures. As shown in Eq. (12. 8p the three-point correlators 
are obtained by combining two open-meson propagators, while for the two-point correlators 
we only need one open-meson propagator. 



C. Doubler modes' effect on the correlation functions 



The remnant doubling degeneracy of staggered fermions leads to contributions of scalar 
states, in addition to pseudoscalar states, in correlation functions with external pseudoscalar 
particles. The scalar contamination yields oscillating terms in the correlation functions (72| . 
In this section, we extend the analysis of Ref. 72|, for two-point correlation functions, to the 
three-point functions introduced in Section IIIBi We conclude that the effect of the doubler 
modes on the three-point functions can be removed at leading order in the lattice spacing 
through appropriate fits of the Euclidean time- dependence. 

The doubling symmetry of the original naive action under the transformation 



where 



(TT, 



^9 = n ^^57m ' 

G = {g : g = (/ii,/i2, 
otherwise 



, /il < /i2 < . . .}, 



9/M 



(2.10) 

(2.11) 

(2.12) 
(2.13) 



generates sixteen equivalent species of quarks, referred to as tastes, that can be reduced to 
four by staggering the quark field fzi. Each element of G is a list of up to four indices, e.g., 
(2), (0,3), and (0,1,2,3) are elements of G, as is the empty set 0. Different g's label different 
doubler modes, or tastes. 
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Consider the general three-point function in momentum space, 

Cr,xr.(tx,t,) = Y,(b{xh5q{x) [g(0)ri6(0)g(0)r26(0)] 6(2/)75g(y)) = 

x,y 

/ ^-)75g(P, t.) [g(0)ri6(0)g(0)r26(0)] t,)75g(fc, ty)) , (2.14) 

where Fi x r2 denotes the Dirac structure of the four-fermion operators in fll.Sp . For 
simphcity of notation, we omit the smearing function from the B meson operator and write 
it as 6(x)75g(x). It would be straightforward (but not particularly instructive) to generalize 
the expressions of Eqs. fl2.14p -f l2.20p to include the smearing function. 

For now the bracketed four-quark operator is left in position space and 6, q are the spatial 
momentum-space bottom and strange/down fermion fields. Because of the doubling sym- 
metry, we can integrate over the central half of the Brillouin zone and sum over the spatial 
doublers 

f.7r/2a 



Cr,xr.(t.,t,) = / 7^^7^^KP^'^9..t,)l^^{v^T^9s^t;) [g(0)Fi6(0)g(0)F26(0)] x 

9a, Ss 



7r/2 

h{k + ng>^,ty)'y5q{k + ng>^,ty)), (2.15) 

where Qs denotes a particular spatial doubler mode. Due to the high momentum that is 
imparted to the heavy quark when Qs 7^ 0, such states are far off-shell and have a negligible 
effect on the correlation function. The taste of the temporal modes can now be considered 
by Fourier transforming the light quarks' temporal component, and then again restricting 
the Brillouin zone and summing over the doublers 

Cr r(t t) - r''" ^^,.o.....o.. .2 16) 

Cr.xr.(t.,t,) - ^2^^3 ^2^^3 j^^^^^ ^2.^ ^^^^e (2.16) 

x^S(p,t,)75 [g'(p,Po) + (-l)*^g'(p,Po + vr/a)] 

X [g(0)Fi6(0)g(0)F26(0)]6(fc,t,)75 [q'{k,ko) + {-lY^q'^k^h + n / a)] ). 
With the momentum space spinors /'^ defined as 



/''(p) = ^z757M?~'(p + vr,) (2.17) 



MG9 

SO that 



9'(P,Po) = /'(p,Po), g'(p,Po + vr/a) = i757o/'°(p,i?o), (2.18) 



the three-point function can be written as 



77r^77rT^AKP.tx)l5 fiP.Q + (-l)^*^^^757o/°(P,tx 

-7r/2a l-^^J l^^J ^ L 

X [g(0)Fi6(0)g(0)F2K0)]5(fc,t,)75 \f{k,ty) + i-lph^,^of\k,ty)] \ ,(2.19) 



9 



where the superscript indicates a temporal taste and no superscript is the null taste at the 
center of the Brillouin zone. 

After Fourier transforming, the bracketed four-quark operator has no restrictions on the 
tastes that contribute to it. However, it must be contracted with the external quark fields 
to form the propagators. Because the asqtad action is used, contractions between tastes of 
different types are suppressed to (9(a^af). The three-point function then takes the form 



/ 7^/21 ^iy, I 

-7r/2a (^TTj-i (iTT)^ \ 

X [(/(O) + z757of (0)) TiKO) (/(O) + ^757o/°(0)) V^h{^)\ b{k,ty)j, 
X [f{k, ty) + (-l)*-2757o/°(fe, ty)] ) , (2.20) 

where higher order terms coming from contractions between quarks of different taste give 
terms of 0{a?a1) that are not considered here. The effects of such terms are comparable 
to NLO terms in staggered chiral perturbation theory and need to be considered at that 
order. They give rise to the "wrong-spin" terms discussed below. According to Eq. fl2.20p . 
the leading-order correlation functions have contributions from both the pseudoscalar and 
the scalar states. The latter ones are known as oscillating states, since the sign of their 
contribution oscillates with time. 

The fit ansatz for our correlators must model both regular and oscillating contributions, 
so that we can remove the latter and extract the physical matrix elements. This is done 
using the form 

JVstates— 1 y^i 

Co''{t,.ty)= V Z^Z. _ ^_^yf.+l)a+fa+l)/3^-E.i.-E,t,^ (2.21) 

a,/3=0 V2E„2E;3 

where the sum is over a finite number of states Ai'states- The time in Eq. fl2.2ip and in the 
discussion on fitting in Sec. |IV] is the number of time slices between the initial state and 
the operator, and thus it is a positive number, unlike the time defined in Fig. [2J The 
oscillations in Euclidean time given by the factor (— l){*^+i)"+(*!'+i)'^ reflect the contribution 
from the scalar states in Eq. f l2.20p . The matrix elements of interest are given by the three- 
point amplitude of the ground state a = /3 = 0, Oqq. Analogously, we incorporate regular 
and oscillating contributions to the description of the two-point correlators by using the 
following functional form in the fits 

-/Vstatcs — 1 

Cps{t)= Yl |^a|'(-l)^*+'^" (e-^"* + e-^"(^-*)) , (2.22) 

where T is the temporal size of the lattice. Three- and two-point functions are fit simulta- 
neously in our analysis, as described in Sec. HVl 



D. Improving the heavy-light four-quark operator 

In addition to the discretization errors in the heavy-quark action, the mixing operator 
also has discretization errors due to the difference in the small-momentum behavior of lattice 
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and continuum heavy quarks. In this section, we describe how the lowest order of operator 
discretization errors are removed in our calculation. We first show that the errors start at 
0{ap) and then discuss how the error at this order can be removed by a "rotation" of the 
heavy- quark field. 

To begin, consider the small ap expansion of the spinor for the Wilson-like fermion 

7o sign X sinh Ea — i'jj sin(pja) + L 



-mia/2 



v/2L(L + sinh Ea) 
Z7 ■ pa 



2 sinhmia 



0((ap) 



-w(x,o) 



(2.23) 



where x labels spin and particle vs. antiparticle, p 

1 

2 



(2/a) sin(pa/2), and for the clover 



action L = 1 + rriQa + ^p'^a'^ — cos poa. The continuum spinor has the expansion 



u 



cont 



ix,p) 



Z7 ■ p 

2m 



0{{ap)' 



w(x,o). 



(2.24) 



The mismatch between the small-momentum terms can be easily removed by "rotating" the 
lattice heavy quark as was done to heavy-light bilinear operators in Ref. j4o|. The light 
lattice spinors for the staggered formulation have the same small-momentum behavior as in 
the continuum up to 0((pa)^, (m^a)^) and need not be matched. 



The analysis of Ref. |40| can be generalized from bilinears to four-fermion operators with 
Dirac structure Fi x We can take one of the terms in the contraction of the lattice 
operator 

Wq)^ Kp'b)\Q^lbq^2b\q{Pq), 6(pfe))lat = 

Nqip'q)Nq{Pq)Nb{p'b)NbiPbMp'q)'^iuTiP'bMPq)'^2U^h\Pb) + (additional coutractlous) , 

(2.25) 

where Nq{p) and Nb{p) are normalization factors for the q and b one-particle states. Following 
the Fermilab interpretation in Ref. 40|], we demand that lattice and continuum amplitudes 
match through 0{ap), 



Z (^n(x,p;)Fie- 
x^(x,Pg)r2e 



-am\ /2 



-am\ /2 



1 - 



n ■ Pb(^ 

2 sinh am\ 

h ■ Pba 
2 sinh am\ 



u{x,0) 
u{x,0) 



u{X,P'q)'^l 

+0{{paf) 



1 - 



^1-Pb 
2mb 



u{x,0) X u{x,Pq)'^2 



1 - 



aZDnQn 

h -Pb 



2m;, 



uix,0) 



(2.26) 



where Qn are dimension-seven lattice operators and -D„ their corresponding coefficients. 
These operators and coefficients are straightforward to identify (m^ must be identified with 
the Fermilab kinetic mass, M2 jz^]). There are two dimension seven operators contributing 
to the matching, Qi = qTij ■ DhqV2h and Q2 = qVihqV^'y ■ Db, which will remove the pa 
discrepancy for appropriate values of the Wilson coefficients and the normalization constant. 
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From the relation above, we find 



Di = D2 



sinh am\ 



2am\ 



and 



Z 



(2.27) 



(2.28) 



Here, mi and m2 are again the Fermilab rest and kinetic masses defined in |40 



However, through 0{ap), adding these operators has the same effect as inducing a "ro- 
tation" of the heavy field 

¥{x) = [1 + adij ■ D]b{x) (2.29) 

where 

di = Di = D2. (2.30) 

This removes 0{Aqcd(i) discretization errors in the operator. The leading errors are then 
((Aqcdc^) ) and 0{asAQCY)a). In the way we have set up the calculation, the open- meson 
propagators, Eq. f l2.9p . include the rotation. 



III. MATCHING OF THE LATTICE MATRIX ELEMENTS 

In order to cancel the scheme and scale dependence of the Wilson coefficients in the 
effective Hamiltonian, we must relate the bare hadronic matrix elements of the lattice oper- 
ators in Eq. (11. 3p to a continuum scheme. We perform that renormalization and matching 
perturbatively at one loop. In the lattice part of this renormalization calculation we use 
mean-field improved lattice perturbation theory [Tsf to improve the convergence of the the- 
ory by resumming the tadpole contributions. 

Already at one loop, even in the continuum, the operators in Eq. (11. 3p mix with each 
other under renormalization. To extract the renormalized value of {Of), we use the following 
matching relation 

{OlY'-^'ifi) =C{[l + a,- Cii(/i,m,,am,)](0?)i-* + a, ■ (12(1^, m,, am,) {OlY^'} 

+ {al,asAQCDa) , (3.1) 

where the renormalization coefficients Qj are the difference between the renormalizations in 
the continuum and on the lattice, Qj = Z^°^^ — Zjj^ . The continuum renormalization scale at 
which we perform the matching is /i. The lattice spacing is a and C is a factor which absorbs 



the lattice field normalization conventions. The values of Zff^^ are listed in Ref. 1761 and a 



detailed description of the calculation of the lattice renormalization coefficients will be given 

in Ref. [t^I • Table [TTl lists the tadpole- improved renormalization coefficients relevant for the 

lattice data analyzed in this paper. For each lattice spacing and b quark mass we show 

the infrared (IR) finite part of the Zjfs as well as the corresponding Qj (in the MS-NDR 

continuum scheme). The Qj are IR finite since the IR divergent contributions to and 

Zff^^ cancel in the difference. All the coefficients Z'f* in Table Hi] are between 0.3 and 1, 
''J ' 

which indicates a sensible behavior of the lattice perturbation series. 

In order to apply the matching relation Eq. (13. ip . we need to choose both a scale n 
and a value for the strong coupling constant a^. For the scale, we use the bottom quark 
mass and, in that way, eliminate higher-order logarithmic contributions that come in powers 
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TABLE II. Values of the finite part of the lattice one- loop renormalization coefficients Z^j^, the 
difference of the continuum and lattice one- loop coefficients Qj needed in (j3.ip for the 0.12 fm and 
0.09 fm lattices, and the coupling used in the matching relation. The continuum (MS-NDR) 
scale used in the matching is /i = mi,. 



a (fm) 


arrib 


r^la,t, finite 
^11 


lylat, finite 
^12 


aMS-NDR 


^MS-NDR 
S12 


as 


0.12 
0.09 


2.1881 
1.7728 


-0.726 
-0.945 


-0.325 
-0.369 


0.1998 
0.3041 


-0.312 
-0.268 


0.32 
0.26 



of log(/i/m 
scheme [70| 

a typical gluon loop momentum in this process and can be calculated using the methods 



For the strong coupling constant, we use the renormalized coupling in the 
evaluated at a scale q*, as in Ref. [78|. The scale q* should be the size of 



outlined in Refs. (70|, l79|. Here, we use q* = 2/a which is close to the calculated value for 
heavy-light currents using the same actions we are employing [t^, S^I- This is justified since 
the contributions coming from the current renormalization are larger than the intrinsic four- 
quark contributions 77|. The values of ay are determined from the static-quark potential 



in a manner similar to that described in Ref. 78|| and are also given in Table [TTl 



IV. FITTING METHOD AND STATISTICAL ERRORS 

The correlation functions are calculated at four different time sources t, and then averaged 
over time sources. For the a ~ 0.12 fm ensembles, to = 0, 16, 32, 48, and for the a ^ 0.09 fm 
ensembles, to = 0, 24, 48, 72. The statistical errors in the data and fits decrease with each 
additional time source by approximately what is expected, suggesting that the correlators 
from different time sources are weakly correlated and statistical power is gained by averaging. 

In order to extract the renormalized matrix elements, we tried two methods for the 
correlator fits. In the first method, we fit the bare correlators and combine the results 
afterwards with the matching coefficients in Sec. Illll to get the renormalized matrix elements. 
In the second method, we first apply the matching coefficients to the correlators for each 
configuration and then perform the fits to obtain the renormalized matrix elements. The 
central values are nearly identical with both methods, but the errors are slightly better and 
the fits more stable with the latter, so for the rest of this article we discuss only the results 
obtained with the second method. 



A. Description of the Fitting Method and Stability Tests 

The heavy-quark in the two-point and three-point correlation functions is always rotated 
at the source as explained in Sec. IIIDi For three-point functions, we smear the heavy quarks 



at the sink using a function based on the quarkonium IS wavefunction [81|, |82|. For two- 
point functions, at the sink we either rotate local heavy quarks, or we smear them with a IS 
wavefunction. Smearing greatly improves the overlap with the ground state. The additional 
rotation at the sink is to ensure that the local-local meson correlator is positive-definite. 
The naive light-quark propagator is always local at source and sink. 

The two-point and three-point correlators used to determine the matrix element (Of) on 
a particular ensemble and for a particular choice of valence-mass rriq are fit simultaneously 
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using the Bayesian fitting approach described in Refs. For the matrix elements on 

the coarse ensembles, we find the smallest errors and greatest stability using three correlators 
(two two-point correlation functions and one three-point correlation function): 



- Cp^ in Eq. fl2.22p with local source and local sink. 

- Cps with local source and IS smeared sink. 

- Cq<i in Eq. f l2.2ip with local source and IS smeared sink. 

For the fine ensembles the best results are obtained with only one two point function and 
one three-point function: 

- Cps with IS smeared source and IS smeared sink. 

- Co<i with local source and IS smeared sink. 

The prior central values function as the initial starting guesses for our fits. Hence we 
choose ground-state values guided by our data to help the fits converge. The prior central 
values for the ground-state masses are obtained from effective mass plots. For the overlap 
factors Zq and Zq*^, where superscripts d and IS" denote factors corresponding to local or 
IS smeared sources/sinks, we examine the amplitude of the B meson propagator with the 
exponential of the ground state removed. We do the same for Oqo, where the Z^^ amplitudes 
are accounted for. The prior widths are taken to be large compared with the statistical error 
of the parameters as reported by the fitter to avoid infiuencing our fit results by our choice 
of priors. For the higher states' overlap factors, the prior width is chosen based on the 
expectation that the overlaps should not be larger than the corresponding ground state 
ones. The energy differences have prior central values and widths that allow them to vary 
from Ai?j+i_j = a(-Ei+i — Ei) k. 0.14 — 0.37, where experimental values [l6[ have been used 
as a guide. We checked that the prior widths for all fitting parameters are large enough so 
they do not infiuence the central value of relevant quantities extracted from the fits. 

The same priors are used for all ensembles, except for the masses of the regular and 
oscillating ground states, Eq and E'^ respectively. These parameters are strongly determined 
by the data, and very different at each lattice spacing, therefore the prior choice must also 
be lattice-spacing dependent. Appendix [C] contains a list of the prior central values and 
widths we use in the calculation. 

Statistical errors are estimated with the bootstrap method. Specifically, for each ensemble 
and valence mass, 500 bootstrap ensembles are constructed from the original ensemble by 
sampling with replacement. A fit is then performed to each ensemble. We find that, as long 
as the bootstrap ensembles are larger than ~ 100, the estimated error is independent of 
bootstrap ensemble size. For fitting methodology checks and plotting purposes in Figs. ISHSl 
statistical errors in the parameters are estimated by the average 68% bootstrap error, which 
is defined as half of the distance between the two points at which 16% of the distribution 
has a higher (lower) value. 

Autocorrelations necessarily exist between correlation functions calculated on different 
configurations within an ensemble and can be minimized by binning the data. The autocor- 
relations are observable only in a few ensemble and valence masses. In many ensembles and 
mass combinations, the noise is large enough that the autocorrelations are not observable. 
We choose a conservative bin size of 4. 
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TABLE III. Number of states and time ranges used for each correlator in the fits for both the 
o ~ 0.12 fm and a ~ 0.09 fm ensembles. For the number of states, the first value indicates the 
number of regular states and the second one the number of oscillating states. The labels between 
parentheses in the first column indicate the type of source/sink in that correlator. 



Correlator 


Number of States 


Time Range 


Cps (local/local) 


3+3 


2-20 


Cps (local/lS) 


3+3 


2-20 


Cq, (local/lS) 


2+2 (2 


- 10) X (2 - 10) 



The number of states included in the sum in Eq. f l2.2ip and the time ranges we use in 
the fits are shown in Table lllli The minimum time slice is fixed to be the same for all the 
correlators in the fit, three-point as well as two-point. However, the maximum time is fixed 



separately for the two- and three-point functions. Following Ref. [83[, the number of states 
are determined by first performing the fit using 1+1 states (1 regular state + 1 oscillatory 
state) starting at large time slices, where the higher energy states no longer contribute 
significantly and a good P^r degree of freedom (d.o.f.) is obtained (~ 1). The fit is then 
performed using one lower time slice as the starting time tmin, and this is repeated, reducing 
tmin until the x^/d.o.f. is no longer reasonable, > 1.5. Then an additional pair of states is 
added to the model function, and the process iterated. Once the timeslice t = 2 can be 
included, that number of states is used in our central-value fitsU 

For the three-point function, we fit using A^'states = -^sTatcfT"^ + -^sTatcs'*™^ = 2 + 2 = 4 and 
timeslices t^, ty G [2, 10] for all ensembles. The two-point functions are fit for t G [2, 20] using 
3 + 3 states. The output of these fits successfully describe the oscillations in the correlation 
functions as can be seen in Fig. [31 which shows the typical behavior in one of the ensembles 
analyzed. 

In order to check that this number of states is sufficient, we add more states and examine 
the stability of the fits. Stability plots over numbers of states for the a ~ 0.09 fm ensembles, 
which illustrate the typical behavior of our fits, are shown in Fig. HI The stability of central 
values and errors is very good for A^statcs > 4 in all cases. 



V. CHIRAL PERTURBATION THEORY 

The light sea- and valence-quark masses that are used in our lattice simulations have 
unphysically large values, with our lightest pion mass ~ 240 MeV. To obtain information 
about the quark-mass dependence of the relevant matrix elements, which allows us to ex- 
trapolate our results to the physical masses, we perform our calculation at six sea x six 
valence quark masses, thus including numerous partially quenched data points. In addition, 
the leading-order taste violations, which arise at 0{a'^as), are included in the theory and 
then removed when the extrapolation is performed using rooted heavy-meson staggered chi- 



ral perturbation theory (rHMSxPT) [63|, |6J, |67| . The xPT expression for {OD, as well as 



for the matrix elements of all the other operators in the AB = 2 effective Hamiltonian was 
first described in Ref. jssi for partially quenched Wilson-type quarks in the framework of 

^ Time slices t = and t = 1 contain unconstrained contamination from higher energy states. At i = all 
states contribute because they have the same exponential weighting, and < = 1 is contaminated by higher 
energy states because the degrees of freedom of staggered fermions spread over two time slices. 
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FIG. 3. Comparison of a fit with correlator data for the correlation function CQi{tx,ty) at fixed 
tx = 2-6 (labeled by color) and as a function of ty, for the a ~ 0.09 fm ensemble with quark masses 
0.0124/0.031 and valence quark mass aniq = 0.031. The lines connect the fit function results for 
integer values of tx,ty coming from the same single fit and evaluated at specific t^, and the dots 
are the average over simulation data. Statistical errors on the simulation data are smaller than the 
plot symbols. The fit results describe very well the oscillation in time shown by the data. 

continuum heavy-meson xPT (HMxPT). With staggered fermions, we also must include the 
effects of taste- violating interactions, using rHMSxPT j67| . 

With the four-quark operators, a careful examination of the Fierz properties shows that 
there are additional operators with both wrong taste and spin, i.e, wrong (Fi, F2). As far as 
we know, this property of local heavy-staggered four-quark operators has not been discussed 
in the literature before. The needed rHMSxPT expressions are derived in Ref. js^. We 
became aware of these contributions after our analysis was nearly complete, so we have not 
included them in the chiral fit functions used here. We do, however, estimate the associated 
systematic error on ^ in our error budget (cf. Sec. I VI Cp . Explicit expressions for the chiral 
fit functions used in this work are given in Appendix IM 

The NLO rHMSxPT in Eq. (lAip and subsequent equations in the appendix can be 
schematically written as 

{B,\Ol\B,) = luljlB, = 

Ms^a [1 + (Wg + Tg + Qg) + L^rrig + L,(2m, + nih) + L^] , (5.1) 

where a, Ly, Ls, and La are low-energy constants (LECs) to be determined from fitting 
the data, the factor of M^^ comes from the HQET normalization of states, and the masses 
rriq, mi, and rrih are the light valence, light sea, and strange sea-quark masses, respectively. 
The light sea quarks are treated as degenerate, and the isospin average is used, i.e., m = 
(m^j + mrf)/2. For staggered quarks the taste-nonsinglet pseudoscalar meson masses are split 

Mfj^p = lj{mi + mj) + a^Ap, (5.2) 
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FIG. 4. fsq sJMBq for the two a ^ 0.09 fm ensembles for anig = 0.031, 0.0031 as a function 
of the number of states A'^jtatcs- The fit results for fBq-\/^Bq Bb^ reach a plateau for A'states > 4. 

where rrii and rrij are the quark masses and the sixteen meson masses are labeled by their 
taste representation, p = P,A, T, V, I. The parameters fi and the A^'s are determined from 



lattice calculations for pions and kaons 87l|. Their values are collected in Table HVl 

The chiral logarithms, Wg, Tq, Qq, stem from wavefunction renormalization, tadpole, and 

sunset diagrams, respectively. The explicit expressions can be found in Appendix El 

In this work, we do not include the effects of the hyperfine splitting A* or the light flavor 

splittings 6k defined in the Appendix. The wrong-spin terms contribute to the tadpole and 

sunset diagrams j86| . 

When extrapolating to the physical point, we set the parameters Ap and (which 
describe discretization effects) and the lattice spacing to zero, and set the sea quark masses 
to their physical values, mi — )■ {niu + md)/2, and nih m^. We then obtain {B^\Of\B^) or 
(5° I Of 1 5°) by setting niq = rrid or nis. Thus, it is an extrapolation to the u- and d-quark 
masses and an interpolation to the s-quark mass. 

An additional consideration is that SU{3) NLO xPT may not be valid for data with 
masses as large as the strange quark's. It would be desirable to include NNLO contributions 
to test the validity of the NLO expression, but the effort needed to calculate the NNLO 
logs is prohibitive. It is reasonable instead to test the chiral expansion by including just 
NNLO analytic contributions. Wherever the quark masses or splittings are large enough for 
such analytic NNLO terms to be significant, the NNLO logarithms should be slowly varying 
and well approximated by the analytic terms. We follow this strategy and supplement the 
NLO rHMSxPT expressions with NNLO analytic terms in our fits, with prior constraints 
estimated based on xPT power counting as explained in the next section. 



A. Parametrization of the Chiral Expression 



Dimensionful quantities are extracted first in units of the lattice spacing and then con- 
verted to their physical values using the ri scale |88|, l89|. This absolute scale is defined 



17 



as r^F(ri) = 1.0, where F{ri) is the force between static quarks. In our chiral fits, all 
parameters are first converted to units of ri by multiplying by the relative scale ri/a. The 
values for ri/a on every MILC ensemble used in our calculation are listed in Table [H After 
the chiral-continuum extrapolation, we convert from ri units with a physical value of ri. 
We take the result obtained by combining the 2009 MILC determination of ri/^r [90[ and 



the PDG value of /^r [16|. Following Ref. [59Lthe error for ri is determined by averaging 



the MILC value with the HPQCD value in [91| and then by inflating the uncertainty to 
take conservatively into account the possible correlations coming from the use of the same 
configurations in both determinations. The final value we use is ri = 0.3117(22) [s^. The 
error associated with ri has a very small effect on the dimensionless quantity C,- 



The dominant lattice artifacts to take into account in our rHMS^PT expressions are 
expected to be taste violating contributions of O(a^a^), since the O(a^as) taste-violating 
effects are absent for asqtad quarks. We parametrize these effects in our fits by defining a 
quantity A^, which is the ratio of the size of taste violations on lattices with spacing a to 

2 

0.12 fm 



those on the a ~ 0.12 fm lattices. Thus A?,,r,r = 1 and 



a2 _ (Q^sQ )o.09 fm ^ op: /p: q\ 

^0.09 fm = / 2 2^ ~ ^■'^^ ■ y^-'^) 

[a^a jo.l2 fm 



The NLO rHMS^PT function in Eq. flS.ip can then be rewritten as 



/3, ^J^{B,\01\B,)/Mb, = Ib^Mb^Bb, 



(5.4) 



where (3^ = ^Ja. The masses Mj^ are defined in Eq. (15. 2p . but here we disregard corrections 
in the masses since they can be absorbed by a redefinition of the low energy constants at 
higher order in the chiral expansion. To the NLO expression above we add, inside the square 
brackets, the allowed NNLO analytic terms, which contribute with seven more unknown 
LECs 

Q^M^^ + Q^{2Ml + Ml^f + Q^Ml^ml + M,\) (5.5) 
+g4(2M,t + M,\) + Pi^X, + P^Alml + Ml,) + P^Al . 

In the above expressions, we have suppressed the factors of ri for simplicity. They can be 
deduced using dimensional considerations. In these expressions, which are the ones we use 
as fit functions, we write the analytic terms for convenience as functions of the pseudoscalar 
masses Mjj rather than the quark masses. 



The ratio ^ can be extracted by first interpolating /3g to vriq = rus and extrapolating to 
= separately according to expressions (15.41) and (15. 5p . and then forming the ratio 
Ps/ Pd- Alternatively, one can consider the ratio of chiral expressions and expand up to 
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TABLE IV. Inputs for the priors of the free parameters and for the fixed parameters in the fits. 
The NLO low energy constants L^, Ls, and La are not constrained in the fits. The parameter s 
is given by the quantity l/(87r^(ri/7r)^). We do not consider errors on the slope firi or the taste 
splittings rfa^Ap because those have negligible effect on the final results. In the right hand side 
table, the two last columns correspond to lattice spacings a ~ 0.12 fm and a ~ 0.09 fm. See the 
text for explanations of the choices of parameters. 



Fit parameters (central valueiwidth) 





a ^ 0.12 fm (a ^ 0.09 fm) 


Input (fixed) parameters 




1 ± 1 


a r 


a 0.12 fm a 


^ 0.09 fm 


gB*B7T 


0.51 ±0.20 




0.2106 


rla?5y 


0.0 ±0.07 (0.0 ± 0.07) X 0.35 




6.234 


6.382 




-0.28 ± 0.06 (-0.28 ± 0.06) x 0.35 










Lv 


unconstrained 




0.439 


0.152 


Ls 


unconstrained 


rfa'^Ar 


0.327 


0.115 


La 


unconstrained 


r\a^AA 


0.205 


0.0706 


Ql-4 


0±S2 

^ , o 


r\a^Ai 


0.537 


0.206 



Q,-w,-r,)±^(M^%,-M;j 

±Qi(M,V - MJJ ± Q^{Ml,^, - Ml^){2Ml + M,\) 
+Pi(M2^,-MjX' (5-6) 

with rHqi fixed to the value closest to m^. In Eq. fl5.6p we disregard the NNLO terms coming 
from squaring the NLO terms in the denominator, since they are not necessary to obtain 
good fits and they are difficult to disentangle from those already included. We can then 
interpolate/extrapolate to m^/ = and = m^^. We call these two strategies for the chiral 
and continuum extrapolation of ^ the indirect and direct methods, respectively. Many of 
the fit parameters cancel in the chiral expression for ^ in Eq. (15. 6p . improving the reliability 
and stability of the fits. In addition, discretization errors of Oias^qpua, (Aqcdo)^) from the 
heavy-quark action that are not included in the chiral perturbation theory, partially cancel 
in the ratio. We thus choose this method as our preferred fitting strategy. 



NNLO to obtain 



B , V^^,' 1 



B. Results from the chiral fits 

In order to perform the chiral fits, we first create 200 bootstrap samples of /3g for each 
sea- and valence-quark mass combination from the two- and three-point correlator fits. The 
bootstrap data is then fit to the chiral expression using Bayesian techniques. The fits are 
simultaneously performed to all ensembles in Table [H 

The input and fit parameters are set as in Table IIVI We do not impose any constraint on 
the NLO low energy constants. For the NNLO LECs we use prior widths based on a simple 
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TABLE V. Results from the rHMSxPT fits. Errors are only statistical and obtained from 200 
bootstrap samples. For a full discussion of systematic errors, see Sec. IVIl 



Ansatz 


xVd.o.f. 


Direct 




f. C, Indirect 


NNLO 
NLO 


0.45 
0.78 


1 9«o+0.035 
i.ZUO_Q 044 

1 90/1+0.018 
i./:«4_ooi6 


0.23 
0.49 


1 9P;p;+0.034 

1 9«9+0.008 
i.ZUZ_o.oi2 



power counting argument. The NLO analytic terms should be of magnitude similar to the 
NLO logs, which are ~ m^/(87r^/^) = s(rim^)^, (with s = l/(87r^(ri/^)^)). Hence, the 
NNLO terms are ~ (m^/(87r^/^))^ = s^(rim^)^. The taste-violating hairpin parameters, 



6y and 5^, were also determined from lattice calculations for pions and kaons in Ref. (62 
We constrain the parameters and 5^ in our fits using the results of Ref. j62[ as prior 
central values and widths. We also take the effective coupling of the B*Btt interaction, 
9b* Bit, as a fit parameter in our analysis. The prior central value and width we use for this 



parameter, shown in Table llVt covers the main ranges of determinations of gB*Bn 92N98 
as discussed in Ref. |57|. A more recent, precise value of gB*BTT, obtained with Nj = 2 + 1 



domain wall fermions and static b quarks [99[, was not yet available when this stage of 



the analysis was carried out. Nevertheless, the result obtained by the authors in Ref. [99 



0.449 ± 0.047 ± 0.019, falls well within the prior central value and width considered here. 



For the pion decay constant, we use the PDG value, /,r = (130.41 ± 0.20) MeV [16 



The fit results for ^ using different ansatzes for the fitting function and the direct and 
indirect methods explained in Sec. IV Al are listed in Table IVl The results and errors obtained 
using the direct and indirect methods agree very well, especially when NNLO terms are 
included. This constitutes a good check of how well our results are encompassing higher- 
order terms in the chiral expansion, which are different in these two methods. 

In Figs. [5] and Ej we show the NLO and NNLO fit results for ^ from the direct method 
as a function of the light valence mass in ri units, rirrig. The top plots in both figures show 
only the full QCD points, rrig = rrii, while the bottom plots show all (partially quenched) 
data included in the fits (see Table [T]). The fit curve is the same in both plots of each figure. 
The black lines show the results of the fit in the continuum limit, after the dominant lattice 
artifacts are removed using rHMSxPT, and after interpolating the physical sea and valence 
strange-quark masses to the physical value, as a function of the valence light-quark mass. 
The black point is our result for at the physical masses, and includes statistical errors. 

From the spread of data in the bottom plots of both figures (same data), one can see 
that the light sea-quark mass dependence is mild; all different sea-quark masses (squares or 
triangles at a particular axis value) agree within one statistical a. The discretization errors 
are also small, as can be seen in both the data and the extrapolation lines in the upper plots. 

We obtain fits that match the data well and have good x^/d.o.f. with only the inclusion of 
NLO terms, as shown in Fig. [51 When we add the NNLO terms, the central values for ^ are 
also within one statistical a, although errors are significantly larger. This is to be expected, 
since the NNLO LECs are poorly known. Related to this is the fact that the x^/d.o.f. for 
the NLO fits are larger than for the NNLO fits. At NNLO, we are including extra degrees 
of freedom with large prior widths that are poorly determined by the fit, so, in practice, we 
are dividing the same by a larger number of degrees of freedom. In fact, the NNLO fits 
seem to give a slightly better description of the data, as can be seen in the full QCD plots. 
The chiral extrapolation for is also milder in the NNLO case. Based on these arguments 
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FIG. 5. Fits results using the NLO rHMSxPT, first line in Eq. ([52]) • The (black) star is the 
physical value of ^ in both plots. The top plot shows only the full QCD data, while the bottom 
one shows all the data included in the fits. The (green) squares and (red) circles and lines in the 
upper plot represent the 0.12 fm and 0.09 fm data and fit results respectively. In the bottom plot, 
each color (symbol) labels a different ensemble. 



and, as mentioned above, the fact that direct and indirect fits agree better at NNLO, we 
choose the direct NNLO fit for our central value and statistical error. The systematic error 
associated with our choice of fit function is discussed in Sec. IVI D[ 
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FIG. 6. Fits results using the NNLO rHMSxPT in Eq. ()5.6p . The (black) star is the physical value 
of ^ in both plots. The top plot shows only the full QCD data, while the bottom one shows all 
the data included in the fits. The (green) squares and (red) circles and lines in the upper plot 
represent the 0.12 fm and 0.09 fm data and fit results, respectively. In the bottom plot, each color 
(symbol) labels a different ensemble. 



VI. ERROR ANALYSIS 



In this section, we discuss all sources of systematic uncertainty affecting our calculation 
of C,. The systematic errors have to be added to the statistical uncertainty listed in Table IVt 
which also encompasses our imperfect knowledge about chiral parameters such as gB* b-k and 
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V,A- 



A. Heavy-quark mass uncertainty 

The mixing parameters depend on the b quark mass used in our simulations tlu'ough 
the hopping parameter which is tuned so that the kinetic meson mass, M2, agrees with 
experiment. The dispersion relation for a heavy particle can be written for low-momentum 

as 



E{p) = Ml + 



PL 

2Mo 



6 



E 



Pi 



VP 



2\2 



(6.1) 



where enters into the definitions of Mi and M2 |40|. and the deviation of M4 from 
M2 capture lattice artifacts. We calculated two-point functions for the pseudoscalar and 
vector mesons at several momenta and extract the energy, E{p), for each particle at each 
momentum. A fit to the dispersion relation then determines M2 and the spin average of the 
results is taken. The Kb value is then adjusted until M2 agrees with the spin-averaged Bs 
meson mass. 

In this work, we use the values for Kb on the a ~ 0.12 fm and a ~ 0.09 fm ensembles 
tuned this way in Ref. 
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The error in the determination translates into a systematic 
error in the mixing matrix elments. However, in the ratio ^ the effect of the uncertainty is 
minimal, since the corrections go in the same direction in both denominator and numerator, 
and, thus, largely cancel. In addition, most of the remaining dependence is encoded in the 
decay constants rather than in the bag parameters, which are very insensitive to the exact 
values of the quark masses. In Ref. 59| , we studied the decay constants with the same choice 
of actions, parameters, and configurations as here. We expect systematic errors to be very 
similar in both analyses. We therefore adopt the error due to the uncertainty in the b quark 
mass obtained in Ref. jsof for the ratio of decay constants fBs/fB^i namely, 0.4%, as a good 
estimate of this systematic error for ^. 



B. Higher-order effects in the perturbative matching 

The most straightforward and conservative way to estimate the effects of the missing 
higher order terms in the perturbative matching is to assume two-loop coefficients of order 
1 and to multiply the central value by a1 = ay {2 /a). This estimate gives an error ~ 5% for 
/eVBe on the a ~ 0.12 fm lattices and ~ 3.6% on the a ~ 0.09 fm lattices, becoming the 
main source of uncertainty for this quantity j33|. If there was no mixing between (Oi) and 
{O2) under renormalization, there would be an exact cancellation of the renormalization 
coefficients for the ratio ^ = (/b^ Bb^) / {fBci\/ Bb^), as long as the valence light-quarks are 
taken to be massless in the renormalization calculation. The mixing under renormalization 
prevents this exact cancellation from happening, but the renormalization corrections in 
the ratio are still largely suppressed, by a factor of {O2) / {OD — with respect 

to those for a single matrix element. We estimate this suppression factor via the ratio 
{nis — md) /Aqcb and multiply the perturbative error for fBy/Bs given above by it. As 
a result, the perturbative matching uncertainty for ^ from this estimate is 0.2-0.5% (for 
Aqcd = 700 MeV). 
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The ratio ^ changes by 0.2% when the one-loop renormahzation is omitted entirely, sup- 
porting our power-counting argument. Another way of estimating 0{a^) effects is by varying 
the scale q* at which ay is evaluated. If we change q* from our central value of 2/ a to 1/a 
and 3/a we find that the extrapolated ^ changes between 0.2-0.4%. 

Since the initial estimate yields the largest uncertainty, 0.5%, we take this as the error 
associated with the missing higher order terms in the perturbative renormahzation. Hence, 
this source of uncertainty is subdominant in our determination of ^. 

C. Mixing with wrong-spin four-fermion operators 

As mentioned in Sec. El there are contributions at NLO in rHMSxPT originating from 
the mixing of (Oi) with the matrix elements of four-fermion operators of different spin and 
taste. We have omitted these contributions from our chiral fits, because we discovered these 
terms after this stage of the analysis was complete. From Figs. |S]and[21 one can see that the 
effect of the wrong-spin mixing is unlikely to be very large, perhaps being mostly absorbed 
into the LECs. 

We cannot include the effects of the wrong-spin contributions, because they require the 
matrix elements of O^, which we have not computed here. Fortunately, however, we have 
started a more comprehensive analysis of B-B mixing on a larger set of higher-statistics 
ensembles, including O3. We have added the wrong-spin operators to that analysis and 
find that their inclusion tends to increase the slope of the continuum extrapolated chiral fit 
function for (Oi) and, hence, ^. For example, taking priors and widths similar to those in 
Table IIVI we find a 2% increase in ^, while for other reasonable choices of the priors the 
variation is not larger than 3.2%. We add a 3.2% systematic error to account for the missing 
terms in our chiral extrapolation functions. 

D. Chiral-extrapolation systematics and light-quark discretization 

The errors due to the choice of fit ansatz and light-quark discretization effects cannot be 
disentangled, because every fit ansatz necessarily treats the discretization errors differently. 
So any estimate of the systematic uncertainty associated with the choice of ansatz also 
accounts for the light-quark discretization errors left over after removing the dominant ones 
using rHMSxPT. 

In Fig. [71 we show the distribution of values for obtained with the NNLO direct fit 
for the 200 bootstrap samples analyzed. We check that 200 bootstrap samples is enough to 
obtain a (nearly) Gaussian distribution, as can be seen in the plot. With the goal of testing 
our choice of functional form and the error associated with the truncation of the chiral 
series, we perform fits with only two of the three NNLO terms, omitting each one in turn. 
All fits give good values of x^/d.o.f. and p value. The values of ^ obtained are scattered 
around the distribution in Fig. [7] but always within 1.5 statistical a of the central value. 
This consistency, together with the fact that the NNLO LECs are not well determined by 
the fit, indicates that the statistical error already accounts for the possibility of having one 
of the unknown constants equal or close to zero. If we inflate and symmetrize our statistical 
errors to ±0.047, we cover the spread of results from the different fitting functions tried 
(including the NLO one). We take this value as our estimate of both statistical and chiral 
systematic uncertainties. 
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FIG. 7. Histogram of the distribution of values of ^ obtained from the 200 bootstrap samples. 
The red line is a Gaussian distribution corresponding ^ = 1.268 ± 0.047 (the NNLO result with 
augmented errors as explained in the text). 



TABLE VI. Input for the physical light-quark masses use d in the chiral extrapolations. These 
values were determined by the MILC collaboration 87, lOol ]. Physical values are found from chiral 



fits that have been extrapolated to the continuum, but masses are still in units of the a ~ 0.09 fm 
lattice spacing. 



Quantity 


Physical 


arus X 10^ 


2.72 (8) 




0.997(35) 


amd X 10^ 


1.40(6) 



An alternative way of estimating the uncertainty in the truncation of the chiral series 
and the fitting function would be taking the difference between the NLO and the NNLO 
fits results. If we add this difference with the statistical errors in Table |V] in quadrature, 
the uncertainty would be slightly smaller than the ±0.047 we are taking as our estimate of 
these two sources of error. 

In the rest of this section, we list the errors associated with the uncertainty of several 
input parameters used in the continuum and chiral fits, that typically can be estimated by 
varying the inputs and redoing the fits. 



1. Light- quark mass uncertainty 



The physical values of the light-quark masses use d for the extrapolations and interpola- 

They are obtained by making 



tions for are determined by the MILC collaboration |87l . ll00 



the charged pions and kaons take on their physical values after removal of electromagnetic 
effects and are listed in Table |VT1 
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The error on ^ due to the hght-quark mass uncertainties is obtained by individually 
varying each quark mass within this uncertainty and repeating the preferred chiral fit and 
extrapolation. The central values arrived at using each mass variation are compared to the 
results of the fit which used the central values for the masses, and their differences are added 
in quadrature. This gives a total systematic error due to the light-quark mass uncertainties 
of 0.5% for e 



2. Uncertainty in the scale ri 

The value of ri used in this analysis to convert from lattice to physical units, as described 
in Sec. IV A[ is ri = 0.3117(22) fm. The results discussed in previous sections are obtained by 
fixing ri to its central value. In order to estimate the uncertainty due to the error in ri we 
change ri by ±0.0022 fm and all parameters that depend on the physical ri are appropriately 
adjusted. The uncertainty in scale gives a systematic error of 0.2%, which is very small due 
to the fact that ^ is a dimensionless quantity and the scale only enters in the normalization 
to lattice units of the chiral corrections (l/(/7rri)^) and indirectly via the tuning of the quark 
masses. 



E. Heavy-quark discretization effects 

The discretization errors associated with our choice of heavy-quark action to simulate 
the bottom quarks can be described in terms of the difference in the lattice and continuum 
Wilson coefficients of higher dimension operators in the HQET expansion. Those come 
from two sources: the mismatch between continuum and lattice in the Lagrangian and the 
mismatch in the four-fermion operators whose matrix elements yield JbVBb and, thus, ^. 
For a particular operator Qi the error can be written in terms of the usual power cou nting 
magnitudes times functions that reflect the particular moa dependence of the action 56, 10l| 



error,/ = Zifi{moa) (aAqcD)^" , (6.2) 

where Sj = dim Qi — 4 for Lagrangian operators Qi of dimension 4 and 5, and Si = dim Qi — 6 
for four-fermion operators Qi of dimension 7 and 8, and the Zi are constants. The functions 
fi{moa) can be deduced from references [i^, 102l | and were discussed in detail in 56|| for the 



form factors parametrizing B — nlu decays and in [59| for heavy-light decay constants. A 
detailed study of these corrections for the matrix elements of all the operators contributing 
to neutral B mixing in the SM and beyond will be presented elsewhere [t^]. Here we only 
summarize the sources of the different corrections for (Oi) and ^. The explicit form of the 
different functions fi{moa) can be found in Appendix [Bl 

From the Lagrangian, there are 0{o?) errors and 0{asa) errors which are identical to 
those in Eqs. (A12) and (A19) in Ref. [59|. They are proportional to the functions /^(moa) 
in Eq. (IBip and /^(moa) in Eq. (]B2p of Appendix [Bl respectively. 

From the four-fermion operators, we have 0{a^) errors coming from higher order cor- 
rections to the rotation relation Eq. (12.291) . They are generated by the mismatch between 
lattice and continuum coefficients of the operators qTD'^b, qTi'S ■ Bh and qVcx ■ Eh in the 
same way as in Eqs. (A13) and (A14) of Ref. [s^, but with an extra overall factor of two due 
to the fact that we have two heavy fields in our leading-order operator, not one. These cor- 
rections are proportional to the functions fxijnoo) and fyiiTioa) in Eq. fIBSp of Appendix [Bj 
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TABLE VII. Heavy-quark discretization errors given in percent. 



V^UlllllUUllUll 


f. 

Ji 




Error /^v^ (%) 
(coarse,fine) 


Error ^ (%) 
(coarse, fine) 


0{a^) Lagrangian 


fE 


2 


(0.28,0.16) 




0{asa) Lagrangian 


fB 


2 


(0.96,0.58) 




O(a^) Operator 


fx 


4 


(1.29,0.74) 






fy 


2 


(0.23,0.18) 




O^aso) Operator 


h 


3 


(1.32,0.75) 




Total error 






(2.1,1.3) 


(0.2,0.1) 



respectively. 

The last contribution, and the least straightforward, comes from the 0{asa) corrections to 
the four-fermion operators. In principle, to subleading order, there are a basis of twelve new 
local operators in the effective Hamiltonian. However, using symmetry constraints, Fierz 
transformations and rewriting some combinations as total derivatives, only five independent 
operators remain Because we separate the temporal and spatial parts of the 



operators through this analysis, the total temporal and spatial components of those five 
operators should be compared with the temporal and spatial parts of the leading operators. 
As explained in Ref. [t^], this produces a difference of 3/3(moa), where /s is given in Eq. (lB4p . 

In Table IVIIt we list all the contributions together with the functions /j, and the pro- 
portionality constant Zi in Eq. (16. 2p . We also list the numerical values of the different 
contributions to the heavy-quark discretization error in percentage for /ba/Bb. In order to 
get the numerical results, we use Aqcd = 700 MeV and = ay(2/a) listed in Table HTl 
The final heavy-quark discretization error for JbVBb is 1.3% for the a ^ 0.09 fm ensembles 
and 2.1% for the 0.12 fm ones. 

These errors largely cancel in ^. The effect of the cancellation on the error can be 
estimated by multiplying the errors in JeVBe by a factor of (m^ — frid) /Aqqb which gives 
a final heavy-quark discretization error for C, of 0.2% for the coarse lattice and 0.1% for the 
fine lattice. This agrees well with the estimate of this type of error for the ratio fs^ / fs |59|. 
~ 0.3%, using a very similar set of data and statistics. The strategy followed in Ref. |59| 
differs from the one described here. In that paper, terms of the form in (16. 2 p were directly 
added to the chiral and continuum extrapolation fitting functions with a coefficient of order 
one to be determined by the fit. Ultimately, we would like to employ that strategy also for 
B^-B^ mixing studies. For this work, however, we simply take the larger estimate from the 
ratio fss/ fs as our estimate of the uncertainty in C, due to heavy-quark discretization errors. 



F. Finite volume corrections 



In order to evaluate the finite volume corrections in our calculation, we follow the pre- 



scription in Refs. [67| and |104| . The MILC lattices are large enough in the time direction 
that it can be treated as infinite to a very good approximation, so we are interested in cor- 
rections due to finite spatial volume only. They are estimated by replacing infinite-volume 
integrals in the chiral expression with finite sums over the spatial momentum. 

Including finite volume corrections in the chiral expressions and redoing the fits reveals 
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TABLE VIII. Complete error budget and total error for the mixing parameter ^. All errors are 
given in percent. 

Source of uncertainty Error (%) 

Statistics © light-quark disc. © chiral extrapolation 3.7 

Mixing with wrong-spin operators 3.2 

Heavy-quark discretization 0.3 

Scale uncertainty (ri) 0.2 

Light- quark masses 0.5 

One-loop matching 0.5 

Tuning Hf, 0.4 

Finite volume 0.1 

Mistuned coarse uq 0.1 

Total Error 5.0 



negligible errors, < 0.1%. 



G. Tuning of the tadpole parameter uq 

The tadpole improvement factor uq is a parameter of the gauge and asqtad staggered 
(sea) quark action and is determined from the fourth root of the average plaquette. The 
tadpole improvement factor also enters into the valence light and heavy quark actions. On 
the a ~ 0.09 fm ensembles, the valence quarks are generated with the same values of Uq as 
the sea. However, on the a ~ 0.12 fm ensembles, the valence quark actions use values of uo 
obtained from the average link in Landau gauge instead. The differences between the values 
of Uq obtained with the two methods is around 3-4%. 

The effect on fBs/fs of the mismatch between uq values in the valence and the sea sectors 
of the a ~ 0.12 fm ensembles was estimated to be < 0.1% in Ref 591. Since this is much 



smaller than the errors due to statistics, chiral fits, and continuum and chiral extrapolation, 
we take this estimate as our error on ^. 



VII. DISCUSSION OF RESULTS AND FUTURE IMPROVEMENTS 

The error budget for the 5'f/(3)-breaking mixing parameter ^ described in the previous 
sections is summarized in Table IVnTl For the first error in the table we prefer not to attempt 
to disentangle the statistical, light-quark discretization, and chiral extrapolation errors since, 
as explained in Sec. IVI D[ the lack of knowledge about the LECs at NNLO makes a reliable 
separation impossible. Our final result is 

^ = 1.268 ±0.063. (7.1) 

The total uncertainty is dominated by the combined statistical, light-quark discretization, 
and chiral extrapolation error, and the uncertainty associated with the wrong-spin operators 
in the chiral-continuum extrapolation. 



28 



Combining our result in Eq. (17. ip with the avera ges of the experimentally measured value s 
of the mass differences AM^ = (0.507±0.004) ps"^ [l6| and AM, = (17.69±0.08) ps-HToH 
and the meson masses M^o = (5366.0 ± 0.9) MeV and M^o = (5279.5 ± 0.5) MeV [3 



we 



quote a value for the ratio of the CKM matrix elements 

Vtd 



ts 



0.216 ±0.011, (7.2) 



assuming no new physics in B^^^yB^^^-^ mixing. The error includes the uncertainties in the 
5- meson masses and mass differences but it is strongly dominated by the error in ^. 

We can also take our result for ^ and combine it with the value of the decay constant 



ratio Jbs/ fsa = 1-229 ± 0.026 calculated by our collaboration [59(| to determine the ratio of 
bag parameters 

^ = ^^(13l\ =1.06 ±0.11. (7.3) 

The two results for and fss/fBa correlated, but the statistical analyses were done 
indepedently so we cannot include the correlations in calculating the uncertainty in the ratio 
of bag parameters. Therefore the error shown in the result of Eq. (17. 3p is overestimated. 
However, as part of the future work, we plan to perform a common analysis of matrix 
elements and decay constants, from which we will be able to account for correlations in 
extracting the value of the bag parameters and thus greatly reduce the error in (17.31) . We 
will do the same for the individual bag parameters corresponding to all the operators in the 
basis in (11.30 . 

Our result for ^ in Eq. (17. ip is in good agreement with the HPQCD value obtained in 



Ref. j32|, ^ = 1.258(33). Note, however, that HPQCD did not estimate the effects of the 
wrong-spin operators that appear in the complete NLO chiral expression, so the full error 
in their result may be somewhat bigger than what was quoted. The agreement of these two 
determinations of ^ provides an excellent check of the methodology and systematic error 
study in both analyses. In addition, it helps to increase the confidence in the robustness 
of lattice results for a parameter of great importance in phenomenological studies. In this 
article, we have established and tested the methodology to apply to broader studies of B^ 
mixing with the same lattice formulations for light and bottom quarks as used here. 

Statistical errors could be reduced significantly by expanding the analysis to include 
the full set of available configurations (approximately 2000) at each of the a ~ 0.12 fm and 
a ~ 0.09 fm ensembles. The current runs of our collaboration on the extended ensembles are 
also implementing sources located at a random spatial and time location to reduce further 
the statistical errors. We expect a reduction of the statistical errors by about a factor of 
two. 

The other dominant error of our calculation, the omission in the rHMS^PT analysis of 
terms generated by wrong-spin operators, will be eliminated when a complete analysis is 



done with the full rHMS^PT expressions [86|. A result for ^ that properly includes the 



wrong-spin terms requires the calculation of the continuum matrix elements not only of the 
operator Oi as we have done in this work, but also of O2 and O3, and simultaneous chiral 
and continuum extrapolations of all three matrix elements. 

The discretization errors, related to both heavy and light quarks, will be reduced in a 
straightforward way by simulations at smaller lattice spacing, i.e., on the a ~ 0.06 fm and 
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a ~ 0.045 fm MILC lattices. The reduction of both statistical and discretization errors will 
also yield cleaner and more accurate continuum and chiral extrapolations. Including data 
at smaller lattice spacings will also reduce the uncertainty associated with the perturbative 
matching from the reduction of = ay(2/a). Although not relevant for the reduction of 
the total error in ^, this will be important in the determination of the matrix elements (Oi) 
themselves. 

Similarly, although the uncertainty associated with heavy-quark discretization effects is 
a subdominant source of error in the determination of ^, it is one of the main errors in the 
determination of (Oi) js^]. In order to have a more reliable, data-driven estimation of these 
effects, in our on-going analyses we plan to employ the strategy used in Ref. j59|, in which 
terms like the ones in (16. 2p are included in the chiral-continuum extrapolation fitting form 
with free parameters to be determined from the fit. 

Our new analysis, which incorporates the improvements mentioned above, includes the 



study of the matrix elements of all five operators that contribute to if^^^^ 29| . This will 
allow not only the precise SM determination of AMg^^, A^s,d, and ^, but will also provide 
the nonperturbative inputs needed to put constraints on BSM models using experimental 
data on mixing and related observables. 
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Appendix A: Staggered Chiral Perturbation Theory for B^-B^ mixing 

In this appendix we describe the functional form we use in the chiral and continuum 
extrapolation of the matrix elements {B^\0\\B^). Further discussion, as well as complete 
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NLO rHMSxPT expressions for {B \Of\B^) with i = 1,...,5 and the corresponding bag 



parameters can be found in [86 



At NLO in rHMSxPT and at first order in the heavy-quark expansion we use 



WlOtlB') = a[l 



>v,5 + m-. 



% + Qc 



+ LyTUg + Ls{2mi + ruh) + LaO^ 



(Al) 



a, Ljj, Ls, and La are constants to be determined from the fits to lattice data. The quantities 
in script for the partially quenched 2+1 (m„ = rrid 7^ nis) case are 



(A2) 



5X 



Xi 



Xv 



(A3) 



p. 
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4fH{MS};W)(^' 



E Dli\{{Mf^}-{,i])n 



(A4) 



In the equations above, the index p runs over the taste representation (P, A, T, V, I) with 
degeneracies A''^ {Np =1,4,6,4,1, respectively), and S runs over the sea flavors u, d, s. The 
meson X is made of two light valence quarks q, and mx is its mass. The functions "H and 
X are the integrals defined in Appendix A of |85|. The subscripts on those functions label 
the flavor and taste of the meson masses at which they are evaluated. 



The superscript in T-L is the second argument for that function as defined in |85|. In 
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addition to the hyperfine splitting A* = Mb* — Mb, it includes a light flavor splitting 5^q 
whenever the light flavor of the vector meson in the loop is different from the external flavor. 

The splitting is 5$q = M^o — M^o = 2Ai/i(ms — m^), where Ai and fi are low energy 
constants. The constant Ai comes from heavy quark effective theory, and /i is deflned in 

The residues functions b!^''^^ and D^"''^' in the expressions above are deflned by jGSf 



The mass combinations appearing as arguments of these functions in the 2+1 partially 
quenched theory are 

{MP} = {m^,mx}, 
{MS^}^{m,,m„m,}, 

{^} = {mL,mH} , (A6) 

where is the meson mass made from // sea quarks, and tuh is the meson mass made from 
hh sea quarks. The tastes of these mesons are indicated explicitly in the equations above. 

Since we are not including the effects of the hyperflne splitting A* or the light flavor 
splittings Sk in this work, the functions "H and I appearing in the wave function, tadpole, 
and sunset contributions simphfy to 



3 ,,0 , /M| 



znl. = -3zX, = = - Y^MJ, In f ^ ) . (AT) 



Appendix B: Functions parametrizing heavy-quark discretization errors 

In this Appendix we collect the functions fi needed in Eq. (16. 2 p to estimate the heavy- 
quark discretization errors affecting our calculation. For details on the ori gin o f these func- 
tions and the effects of higher- dimension operators in the lagrangian, see 102l |. For further 



details on the application to the estimation of heavy-quark discretization errors in — 
mixing, see (77j . 



0{a?) errors from the Lagrangian. 

(1 + rriQa) — 1 



fE{moa) = ^ 



(Bl) 



moa(2 + moa)(l + mofl) 4(l + moa)" 

0{asa'^) errors from the Lagrangian 

IbM = J (B2) 
2(1 -|- moa) 
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O(a^) errors from the four-fermion operator 



fx{moa) = - 



/y(moa) 



mo a 



2(l + moa) \2(2 + moa)(l + moa) 
2 + inioa + {niQay 



4(1 + moa)2(2 + moa)^ 



(B3) 



O(asa^) errors from the four-fermion operator 



2(2 + moa) 



(B4) 



Appendix C: Prior central values and widths for the correlator fits 

In the table below we collect the prior central values and widths used in the correlator 
fits described in Sec. IIVI The amplitude parameters are defined in Eqs. f l2.2ip and fl2.22p . 
and the energy differences are defined as A_Ej+i j = a(-Ei+i — Ei). 



TABLE IX. The priors with index refer to the ground state. Superscripts d and IS" refer to the 
local and IS smeared sources respectively. Higher energy state priors have indices i and j. The 
prime in E'^ refers to an opposite parity (oscillating) state. 



Prior central value Prior width 



^0 


2.2 


0.5 




0.01 


0.5 


zt 


0.45 


0.45 


zf 


0.01 


1 


Ooo 


0.01 


0.02 




0.01 


0.1 


Eq (0.12 fm) 


1.95 


0.15 


E'q (0.12 fm) 


2.25 


0.15 


Eq (0.09 fm) 


1.65 


0.15 


E'q (0.09 fm) 


1.85 


0.15 


log AEi+i^i 


-1.5 


0.5 




-1.5 


0.5 
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